function F = cal_firm_vars(param, const, eqm)
 
    % Define function
    firmvar = @(var) var(eqm.q_index)';

    
%% Quality, Productivity

    % Quality level
    firm.qlevel = firmvar(const.Qgrid);
    
    % Productivity
    firm.lnz = full(sum(const.lnzQ.*eqm.q_indicator, 2));

    
%% Import probability of exporters/non-exporters
    
    prob_exp_imp        = eqm.prob_exp_imp;
    prob_exp_nonimp     = eqm.prob_exp_nonimp;    
    prob_nonexp_imp     = eqm.prob_nonexp_imp;
    prob_nonexp_nonimp  = eqm.prob_nonexp_nonimp;
    
        
%% Export and import probabilities/intensity        
        
    % export probability
    firm.prob_exp = full(sum(eqm.ExportProbQ.*eqm.q_indicator, 2));    
    firm.prob_nonexp = full(sum((1-eqm.ExportProbQ).*eqm.q_indicator, 2));

    % import probability
    firm.prob_imp = firm.prob_exp.*prob_exp_imp + firm.prob_nonexp.*prob_nonexp_imp;
    firm.prob_nonimp = firm.prob_exp.*prob_exp_nonimp + firm.prob_nonexp.*prob_nonexp_nonimp;
    
    % export intensity
    firm.exp_intst_exp = firmvar(1-eqm.rv_rev);

    % import intensity
    firm.imp_intst_imp = firmvar(1-eqm.rm_spend);
    
    
    
%% Sales

    % Log sales (E, I)
    firm.lnsales_exp_imp      = (param.sigma-1)*const.gamma*firm.lnz + log(firmvar(eqm.Pi11));
    firm.lnsales_exp_nonimp   = (param.sigma-1)*const.gamma*firm.lnz + log(firmvar(eqm.Pi10));
    firm.lnsales_nonexp_imp   = (param.sigma-1)*const.gamma*firm.lnz + log(firmvar(eqm.Pi01));
    firm.lnsales_nonexp_nonimp= (param.sigma-1)*const.gamma*firm.lnz + log(firmvar(eqm.Pi00));

    % Log expected sales (E)
    firm.lnsales_exp = log(exp(firm.lnsales_exp_imp).*prob_exp_imp + exp(firm.lnsales_exp_nonimp).*prob_exp_nonimp);
    firm.lnsales_nonexp = log(exp(firm.lnsales_nonexp_imp).*prob_nonexp_imp + exp(firm.lnsales_nonexp_nonimp).*prob_nonexp_nonimp);


%% Indegree

    % Number of suppliers (E, I)
    firm.indegree_exp_imp       = exp(firm.lnsales_exp_imp).^(1/param.beta_m).*firmvar(eqm.rm_ads.*eqm.C_m1.*eqm.theta_m);
    firm.indegree_nonexp_imp    = exp(firm.lnsales_nonexp_imp).^(1/param.beta_m).*firmvar(eqm.rm_ads.*eqm.C_m1.*eqm.theta_m);
    firm.indegree_exp_nonimp    = exp(firm.lnsales_exp_nonimp).^(1/param.beta_m).*firmvar(1.*const.C_m0.*eqm.theta_m);
    firm.indegree_nonexp_nonimp = exp(firm.lnsales_nonexp_nonimp).^(1/param.beta_m).*firmvar(1.*const.C_m0.*eqm.theta_m);

    % Number of suppliers (E)
    firm.indegree_exp = firm.indegree_exp_imp.*prob_exp_imp + firm.indegree_exp_nonimp.*prob_exp_nonimp;
    firm.indegree_nonexp = firm.indegree_nonexp_imp.*prob_nonexp_imp + firm.indegree_nonexp_nonimp.*prob_nonexp_nonimp;
       

%% Outdegree

    % Number of customers (E, I)
    tempsum = sum(const.phi_v.*kron(ones(1,param.Q),eqm.theta_v'), 1);
    firm.outdegree_exp_imp       = exp(firm.lnsales_exp_imp).^(1/param.beta_m).*firmvar(eqm.rv_ads.*eqm.C_v1.*tempsum);
    firm.outdegree_exp_nonimp    = exp(firm.lnsales_exp_nonimp).^(1/param.beta_m).*firmvar(eqm.rv_ads.*eqm.C_v1.*tempsum);
    firm.outdegree_nonexp_imp    = exp(firm.lnsales_nonexp_imp).^(1/param.beta_m).*firmvar(1.*const.C_v0.*tempsum);
    firm.outdegree_nonexp_nonimp = exp(firm.lnsales_nonexp_nonimp).^(1/param.beta_m).*firmvar(1.*const.C_v0.*tempsum);    

    % Number of customers (E)
    firm.outdegree_exp = firm.outdegree_exp_imp.*prob_exp_imp + firm.outdegree_exp_nonimp.*prob_exp_nonimp;
    firm.outdegree_nonexp = firm.outdegree_nonexp_imp.*prob_nonexp_imp + firm.outdegree_nonexp_nonimp.*prob_nonexp_nonimp;
    
    
%% Wage, Employment

    % Average wage per worker
    firm.wage = firmvar(param.w.*const.wage_grid);    

%{    
    % Nnumber of workers (E, I)
    frac = (1-param.alpha_m-param.alpha_s)*(param.sigma-1)/param.sigma;
    firm.employment_exp_imp         = frac.*exp(firm.lnsales_exp_imp)./firm.wage;
    firm.employment_exp_nonimp      = frac.*exp(firm.lnsales_exp_nonimp)./firm.wage;      
    firm.employment_nonexp_imp      = frac.*exp(firm.lnsales_nonexp_imp)./firm.wage; 
    firm.employment_nonexp_nonimp   = frac.*exp(firm.lnsales_nonexp_nonimp)./firm.wage;

    % Number of workers (E)
    firm.employment_exp = firm.employment_exp_imp.*prob_exp_imp + firm.employment_exp_nonimp.*prob_exp_nonimp;
    firm.employment_nonexp = firm.employment_nonexp_imp.*prob_nonexp_imp + firm.employment_nonexp_nonimp.*prob_nonexp_nonimp;
%}    
    
%% Average log wage of suppliers

    % Unweighted average supplier wage
    avg_unwgt_lnwageS_q = 1./eqm.V.*sum(const.phi_v.*kron(ones(param.Q,1), ...
        (log(param.w)+const.lnwage_grid).*eqm.Vbar),2)';
    firm.avg_unwgt_lnwageS = firmvar(avg_unwgt_lnwageS_q);
    
    % Expenditure weighted average supplier wage
    avg_wgt_lnwageS_q = eqm.theta_m.*eqm.c_H.^(param.sigma-1)./eqm.V.*sum(const.phi_y...
        .*const.phi_v.*kron(ones(param.Q,1), (log(param.w)+const.lnwage_grid).*eqm.Psig),2)';    
    firm.avg_wgt_lnwageS = firmvar(avg_wgt_lnwageS_q);
      
    
%% Fraction of higher quality match
%{
    % Fraction of higher (or same) quality supplier
    flag_S_higher_q = zeros(param.Q, param.Q); %row is buyer, column is supplier
    for i = 1:param.Q
        flag_S_higher_q(i,:) = [zeros(1,i-1), ones(1,param.Q-i+1)];
    end
    frac_S_higher_q = sum(flag_S_higher_q.*const.phi_v.*repmat(eqm.Vbar,param.Q,1),2)'./eqm.V;
    firm.frac_S_higher_q = firmvar(frac_S_higher_q);
    
    % Fraction of higher quality customer
    flag_B_higher_q = zeros(param.Q, param.Q); %row is buyer, column is seller
    for i = 1:param.Q
        flag_B_higher_q(:,i) = [zeros(i-1,1); ones(param.Q-i+1,1)];
    end    
    tempmatrix = const.phi_v.*repmat(eqm.theta_v',1,param.Q);
    frac_B_higher_q = sum(flag_B_higher_q.*tempmatrix,1)./sum(tempmatrix,1);        
    firm.frac_B_higher_q = firmvar(frac_B_higher_q);
%}            

%% Return

    F = firm; 

    
end